#library(diseasemapping)
load("referencepop.RData")
source("formatCases.r")
source("formatPopulation.data.frame.r")
source("getSMR.data.frame.r")
source("getBreaks.r")
source("getRates.r")

#load census data
load("census_data.RData")
ontPop5 <- NULL
ontPop5[[1]] <- cbind(ontPop[[1]][,1:4], 5 * ontPop[[1]][,-(1:4)])
ontPop5[[2]] <- cbind(ontPop[[2]][,1:4], 5 * ontPop[[2]][,-(1:4)])
ontPop5[[3]] <- cbind(ontPop[[3]][,1:5], 5 * ontPop[[3]][,-(1:5)])
ontPop5[[4]] <- cbind(ontPop[[4]][,1:5], 5 * ontPop[[4]][,-(1:5)])
ontPop5[[5]] <- cbind(ontPop[[5]][,1:5], 5 * ontPop[[5]][,-(1:5)])
canadaPopProp1M <- data.frame(t(1000000 * referencepop$POPULATION))
names(canadaPopProp1M) <- names(ontPop[[5]])[-c(1:7,26)]
canadaPopProp1M <- data.frame(year = 1991, canadaPopProp1M)

#load case data
obstotCases <- read.table("TotCases.txt", header = TRUE)
names(obstotCases) <- gsub("Frequency.1", "cumFrequency", names(obstotCases))
names(obstotCases) <- gsub("Percent.1", "cumPercent", names(obstotCases))
obsagesexCases07 <- read.table("cases07.txt", header = TRUE)
obsagesexCases07$age <- 5 * obsagesexCases07$agegrp - 1
obscdCases <- read.table("stuffrevised.txt", sep = "\t", header = FALSE)
colnames(obscdCases) = c("id", "year", "sex", "sgc", "pc", "dco", "site")
obscdCases$county = substr(obscdCases$sgc, 1, 2)
obscdCases$county = factor(obscdCases$county, levels = 37:39, labels = c("Essex", "Lambton", "Middlesex"))

#---
#Ontario and Canada Population
#---

#incidence rates by age-sex group for Ontario population
agesexModel <- getRates(obsagesexCases07, ontPop5[[5]], ~age:sex - 1, breaks = seq(5, 75, by = 10))
agesexRates <- exp(agesexModel$coef)
agesexRates <- agesexRates[order(names(agesexRates))]

#plot incidence rates for Ontario
pdf("age_sex_rate.pdf", width = 6, height = 4, pointsize = 12)
agesexRatesF <- agesexRates[seq(1, length(agesexRates), by = 2)]
agesexRatesM <- agesexRates[seq(2, length(agesexRates), by = 2)]
plot(1:length(agesexRatesM), 1e6 * agesexRatesM, type = "l", axes = FALSE, ylim = c(0,150),
	main = "Age-Sex Incidence Rates in Ontario", 
	xlab = "Age Group", ylab = expression(paste("Rate (", 10^-6, " Cases per Year)", sep = "")),
	cex.lab = 0.8, lwd = 1)
lines(1:length(agesexRatesF), 1e6 * agesexRatesF, col = 1, lty = 2, lwd = 1)
box(lwd = 1)
axis(1, at = 1:length(agesexRatesM), labels = c(paste(seq(15, 65, by = 10), "-", seq(24, 74, by = 10), sep = ""), "75+"),
	tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 150, by = 25), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
legend("topleft", legend = c("Male", "Female"), lty = 1:2, col = 1, cex = 0.8, lwd = 1)
dev.off()

#observed total number of cases for Ontario between 1983 and 2007
year <- 1983:2007
obstotCases83_07 <- obstotCases[obstotCases$diagccyy >= year[1] & obstotCases$diagccyy <= year[length(year)],2]

#expected total number of cases for Ontario between 1983 and 2007
exptotCases86 <- getSMR(ontPop[[1]], agesexModel, regionCode = "year")$expected
exptotCases91 <- getSMR(ontPop[[2]], agesexModel, regionCode = "year")$expected
exptotCases96 <- getSMR(ontPop[[3]], agesexModel, regionCode = "year")$expected
exptotCases01 <- getSMR(ontPop[[4]], agesexModel, regionCode = "year")$expected
exptotCases06 <- getSMR(ontPop[[5]], agesexModel, regionCode = "year")$expected
exptotCases83_07 <- rep(c(exptotCases86,exptotCases91,exptotCases96,exptotCases01,exptotCases06), each = 5)

#plot observed and expected total number of cases for Ontario
pdf("obs_exp_cases_Ontario.pdf", width = 6, height = 4, pointsize = 12)
plot(year, obstotCases83_07, type = "l", lty = 2, axes = FALSE, xlim = c(1983,2007), ylim = c(40,180),
	main = "Total Number of Cases in Ontario",
	xlab = "Year", ylab = "Cases", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(40, 180, by = 20), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
lines(year, exptotCases83_07, lty = 1)
legend("topleft", legend = c("Observed", "Expected"), lty = c(2,1), col = 1, cex = 0.8, lwd = 1)
dev.off()

#standardized mortality ratio for total number of cases for Ontario between 1983 and 2007
smr <- obstotCases83_07/exptotCases83_07

#plot standardized mortality ratio for total number of cases for Ontario
pdf("smr_Ontario.pdf", width = 6, height = 4, pointsize = 12)
plot(year, smr, type = "l", axes = FALSE, xlim = c(1983,2007), ylim = c(0.5,1.1), 
	main = "Standardized Mortality Ratio in Ontario",
	xlab = "Year", ylab = "Ratio", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0.5, 1.1, by = 0.1), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

#expected total number of cases in Canada
stdRate <- getSMR(canadaPopProp1M, agesexModel, regionCode = "year")$expected

#plot expected rate of cases for Canada
pdf("exp_rate_Canada.pdf", width = 6, height = 4, pointsize = 12)
plot(year, stdRate * smr, type = "l", axes = FALSE, xlim = c(1983,2007), ylim = c(5,12), 
	main = "Expected Rate of Cases in Canada",
	xlab = "Year", ylab = expression(paste("Rate (Cases per ", 10^6, " Persons)", sep = "")), 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(5, 12, by = 1), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

#---
#Census Division
#---

#observed number of cases for Essex County, Lambton County and Middlesex County between 1983 and 2007
obsCases <- table(obscdCases$year, obscdCases$county)
obsCases83_07 <- obsCases[rownames(obsCases) >= year[1] & rownames(obsCases) <= year[length(year)],]

#expected number of cases for Essex County, Lambton County and Middlesex County between 1983 and 2007
expCases86 <- getSMR(cdPop[[1]], agesexModel, regionCode = "cduid")$expected
expCases91 <- getSMR(cdPop[[2]], agesexModel, regionCode = "cduid")$expected
expCases96 <- getSMR(cdPop[[3]], agesexModel, regionCode = "cduid")$expected
expCases01 <- getSMR(cdPop[[4]], agesexModel, regionCode = "cduid")$expected
expCases06 <- getSMR(cdPop[[5]], agesexModel, regionCode = "cduid")$expected
expCasesEssex83_07 <- rep(c(expCases86[1],expCases91[1],expCases96[1],expCases01[1],expCases06[1]), each = 5) * smr
expCasesLambton83_07 <- rep(c(expCases86[2],expCases91[2],expCases96[2],expCases01[2],expCases06[2]), each = 5) * smr
expCasesMiddlesex83_07 <- rep(c(expCases86[3],expCases91[3],expCases96[3],expCases01[3],expCases06[3]), each = 5) * smr

#plot observed and expected number of cases for Essex County
pdf("obs_exp_cases_Essex.pdf", width = 6, height = 4, pointsize = 12)
plot(year, obsCases83_07[,1], type = "l", lty = 2, axes = FALSE, xlim = c(1983,2007), ylim = c(0,10),
	main = "Number of Cases in Essex County",
	xlab = "Year", ylab = "Cases", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 10, by = 2), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
lines(year, expCasesEssex83_07, lty = 1)
legend("topleft", legend = c("Observed", "Expected"), lty = c(2,1), col = 1, cex = 0.8, lwd = 1)
dev.off()

#plot observed and expected number of cases for Lambton County
pdf("obs_exp_cases_Lambton.pdf", width = 6, height = 4, pointsize = 12)
plot(year, obsCases83_07[,2], type = "l", lty = 2, axes = FALSE, xlim = c(1983,2007), ylim = c(0,16),
	main = "Number of Cases in Lambton County",
	xlab = "Year", ylab = "Cases", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 16, by = 2), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
lines(year, expCasesLambton83_07, lty = 1)
legend("topleft", legend = c("Observed", "Expected"), lty = c(2,1), col = 1, cex = 0.8, lwd = 1)
dev.off()

#plot observed and expected number of cases for Middlesex County
pdf("obs_exp_cases_Middlesex.pdf", width = 6, height = 4, pointsize = 12)
plot(year, obsCases83_07[,3], type = "l", lty = 2, axes = FALSE, xlim = c(1983,2007), ylim = c(0,7),
	main = "Number of Cases in Middlesex County",
	xlab = "Year", ylab = "Cases", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 7, by = 1), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
lines(year, expCasesMiddlesex83_07, lty = 1)
legend("topleft", legend = c("Observed", "Expected"), lty = c(2,1), col = 1, cex = 0.8, lwd = 1)
dev.off()

#re-organize data for Essex County, Lambton County and Middlesex County
obsexpCases83_07 <- data.frame(year = rep(year, 3), 
	cdname = rep(c("Essex","Lambton","Middlesex"), each = length(year)),
	observed = c(obsCases83_07[,1],obsCases83_07[,2],obsCases83_07[,3]),
	expected = c(expCasesEssex83_07,expCasesLambton83_07,expCasesMiddlesex83_07))
obsexpCases83_07$logexpected <- log(obsexpCases83_07$expected)

#GAM for Essex County
library(mgcv)
gamEssex <- gam(observed ~ s(year) + offset(logexpected), family = "poisson", 
	data = obsexpCases83_07, subset = cdname == "Essex")
summary(gamEssex)

#plot year effect from GAM for Essex County
pdf("gam_year_Essex.pdf", width = 6, height = 4, pointsize = 12)
plot(gamEssex, se = TRUE,
	main = "Year Effect in Essex County", 
	xlab = "Year", ylab = "s(Year,1) (Linear Predictor Scale)", 
	cex.lab = 0.8, cex.axis = 0.8, lwd = 1)
dev.off()

#plot predictions from GAM for Essex County
gampredictEssex <- predict(gamEssex, se.fit = TRUE)
gamlpEssex <- gampredictEssex$fit
gamlplowerEssex <- gamlpEssex - 1.96 * gampredictEssex$se.fit
gamlpupperEssex <- gamlpEssex + 1.96 * gampredictEssex$se.fit
pdf("gam_predict_lp_Essex.pdf", width = 6, height = 4, pointsize = 12)
matplot(year, cbind(gamlpEssex,gamlplowerEssex,gamlpupperEssex), type = "l", lty = c(1,2,2), col = 1, 
	axes = FALSE, xlim = c(1983,2007), ylim = c(0,2.5),
	main = "Predicted Number of Cases in Essex County",
	xlab = "Year", ylab = "Cases (Linear Predictor Scale)", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 2.5, by = 0.5), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

gamfitEssex <- exp(gamlpEssex)
gamfitlowerEssex <- exp(gamlplowerEssex)
gamfitupperEssex <- exp(gamlpupperEssex)
pdf("gam_predict_fit_Essex.pdf", width = 6, height = 4, pointsize = 12)
matplot(year, cbind(gamfitEssex,gamfitlowerEssex,gamfitupperEssex), type = "l", lty = c(1,2,2), col = 1, 
	axes = FALSE, xlim = c(1983,2007), ylim = c(0,10),
	main = "Predicted Number of Cases in Essex County",
	xlab = "Year", ylab = "Cases (Transform Scale)", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 10, by = 2), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

#GAM for Lambton County
gamLambton <- gam(observed ~ s(year) + offset(logexpected), family = "poisson", 
	data = obsexpCases83_07, subset = cdname == "Lambton")
summary(gamLambton)

#plot year effect from GAM for Lambton County
pdf("gam_year_Lambton.pdf", width = 6, height = 4, pointsize = 12)
plot(gamLambton, se = TRUE,
	main = "Year Effect in Lambton County", 
	xlab = "Year", ylab = "s(Year,1.5) (Linear Predictor Scale)", 
	cex.lab = 0.8, cex.axis = 0.8, lwd = 1)
dev.off()

#plot predictions from GAM for Lambton County
gampredictLambton <- predict(gamLambton, se.fit = TRUE)
gamlpLambton <- gampredictLambton$fit
gamlplowerLambton <- gamlpLambton - 1.96 * gampredictLambton$se.fit
gamlpupperLambton <- gamlpLambton + 1.96 * gampredictLambton$se.fit
pdf("gam_predict_lp_Lambton.pdf", width = 6, height = 4, pointsize = 12)
matplot(year, cbind(gamlpLambton,gamlplowerLambton,gamlpupperLambton), type = "l", lty = c(1,2,2), col = 1, 
	axes = FALSE, xlim = c(1983,2007), ylim = c(0.5,3),
	main = "Predicted Number of Cases in Lambton County",
	xlab = "Year", ylab = "Cases (Linear Predictor Scale)", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0.5, 3, by = 0.5), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

gamfitLambton <- exp(gamlpLambton)
gamfitlowerLambton <- exp(gamlplowerLambton)
gamfitupperLambton <- exp(gamlpupperLambton)
pdf("gam_predict_fit_Lambton.pdf", width = 6, height = 4, pointsize = 12)
matplot(year, cbind(gamfitLambton,gamfitlowerLambton,gamfitupperLambton), type = "l", lty = c(1,2,2), col = 1, 
	axes = FALSE, xlim = c(1983,2007), ylim = c(0,18),
	main = "Predicted Number of Cases in Lambton County",
	xlab = "Year", ylab = "Cases (Transform Scale)", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 18, by = 2), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

#GAM for Middlesex County
gamMiddlesex <- gam(observed ~ s(year) + offset(logexpected), family = "poisson", 
	data = obsexpCases83_07, subset = cdname == "Middlesex")
summary(gamMiddlesex)

#plot year effect from GAM for Middlesex County
pdf("gam_year_Middlesex.pdf", width = 6, height = 4, pointsize = 12)
plot(gamMiddlesex, se = TRUE,
	main = "Year Effect in Middlesex County", 
	xlab = "Year", ylab = "s(Year,2.42) (Linear Predictor Scale)", 
	cex.lab = 0.8, cex.axis = 0.8, lwd = 1)
dev.off()

#plot predictions from GAM for Middlesex County
gampredictMiddlesex <- predict(gamMiddlesex, se.fit = TRUE)
gamlpMiddlesex <- gampredictMiddlesex$fit
gamlplowerMiddlesex <- gamlpMiddlesex - 1.96 * gampredictMiddlesex$se.fit
gamlpupperMiddlesex <- gamlpMiddlesex + 1.96 * gampredictMiddlesex$se.fit
pdf("gam_predict_lp_Middlesex.pdf", width = 6, height = 4, pointsize = 12)
matplot(year, cbind(gamlpMiddlesex,gamlplowerMiddlesex,gamlpupperMiddlesex), type = "l", lty = c(1,2,2), col = 1,
	axes = FALSE, xlim = c(1983,2007), ylim = c(0,2.5),
	main = "Predicted Number of Cases in Middlesex County",
	xlab = "Year", ylab = "Cases (Linear Predictor Scale)", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 2.5, by = 0.5), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

gamfitMiddlesex <- exp(gamlpMiddlesex)
gamfitlowerMiddlesex <- exp(gamlplowerMiddlesex)
gamfitupperMiddlesex <- exp(gamlpupperMiddlesex)
pdf("gam_predict_fit_Middlesex.pdf", width = 6, height = 4, pointsize = 12)
matplot(year, cbind(gamfitMiddlesex,gamfitlowerMiddlesex,gamfitupperMiddlesex), type = "l", lty = c(1,2,2), col = 1,
	axes = FALSE, xlim = c(1983,2007), ylim = c(0,10), 
	main = "Predicted Number of Cases in Middlesex County",
	xlab = "Year", ylab = "Cases (Transform Scale)", 
	cex.lab = 0.8, lwd = 1)
box(lwd = 1)
axis(1, at = seq(1985, 2005, by = 5), labels = TRUE, tick = TRUE, las = 1, cex.axis = 0.8)
axis(2, at = seq(0, 10, by = 2), labels = TRUE, tick = TRUE, las = 2, cex.axis = 0.8)
dev.off()

save.image("meso_analysis.RData")
